Banded Dotterel Moult Study

Exploration of dataset

Authors
Affiliations

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Bashar Jarayseh

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Ailsa Howard

South Bay Banded Dotterel Project, Kaikoura, New Zealand

Tony Habraken

Port Waikato Banded Dotterel Project, Port Waikato, New Zealand

Emma Williams

Department of Conservation, Christchurch, New Zealand

Colin O`Donnell

Department of Conservation, Christchurch, New Zealand

Bart Kempenaers

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Published

November 13, 2024

Code
knitr::opts_chunk$set(cache = TRUE)

Explore Bashar’s Dataset (modified 7-Nov-2024)

This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura

Plot the sampling distribution of Ailsa’s photos across the seasons

Code
# 2021-2022 season data
ggplot(data = dat %>% filter(season == 2021 & sex == "F") %>% mutate(score = as.numeric(score), 
                                                                     rings_comb_ = reorder(rings_comb, n_photos, decreasing = TRUE))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ rings_comb_, ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +
  xlab("week") +
  ggtitle("Females from the 2021-2022 breeding season in Kaikoura")

Code
ggplot(data = dat %>% filter(season == 2021 & sex == "M") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +
  xlab("week") +
  ggtitle("Males from the 2021-2022 breeding season in Kaikoura")

Code
# 2022-2023 season data (not as good coverage as the previous season)
ggplot(data = dat %>% filter(season == 2022 & sex == "F") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +
  xlab("week") +
  ggtitle("Females from the 2022-2023 breeding season in Kaikoura")

Code
ggplot(data = dat %>% filter(season == 2022 & sex == "M") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +
  xlab("week") +
  ggtitle("Males from the 2022-2023 breeding season in Kaikoura")

Code
# 2023-2024 season data 
ggplot(data = dat %>% filter(season == 2023 & sex == "F") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2023-07-01"), as.Date("2024-05-01"))) +
  xlab("week") +
  ggtitle("Females from the 2023-2024 breeding season in Kaikoura")

Code
ggplot(data = dat %>% filter(season == 2023 & sex == "M") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2023-07-05"), as.Date("2024-05-01"))) +
  xlab("week") +
  ggtitle("Males from the 2023-2024 breeding season in Kaikoura")

sex-specific differences in breeding plumage

Code
# density plot of sex-specific breeding plumage score distributions
ind_breeding_scores %>% 
  group_by(rings_comb, sex) %>% 
  summarise(max_breeding_score = max(max_breeding_score)) %>% 
  ungroup() %>% 
  group_by(sex, max_breeding_score) %>% 
  summarise(n = n()) %>% 
  bind_rows(., data.frame(sex = c("F", "M"), max_breeding_score = c(7, 3), n = c(0, 0))) %>%
  ungroup() %>% 
  ggplot() +
  geom_bar(aes(x = max_breeding_score, y = n, fill = sex), 
           alpha = 0.5, stat = "identity", position = "dodge", width = 0.5) +
  luke_theme +
  theme(legend.position = c(0.15, 0.9)) +
  xlab("Maximum individual rufous band score during breeding season") +
  ylab("Frequency") +
  scale_colour_brewer(palette = "Dark2", direction = -1,
                      name = "Sex",
                      labels = c("Female (N = 49)", "Male (N = 35)")) +
  scale_fill_brewer(palette = "Dark2", direction = -1,
                    name = "Sex",
                    labels = c("Female (N = 49)", "Male (N = 35)"))

Code
# draw gt table
mod_max_score_table <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = c(estimate),
             rows = 1:10,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = c(estimate),
             rows = 11:13,
             decimals = 0,
             use_seps = FALSE) %>% 
  sub_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = c(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

mod_max_score_table
Banded Dotterel rufous band score Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Male (Sex) 1.01 [0.74, 1.28]
Partitioned 𝑹²
Total Marginal 𝑹² 0.29 [0.17, 0.4]
Sex 0.29 [0.17, 0.4]
Total Conditional 𝑹² 0.51 [0.35, 0.66]
Random effects 𝜎²
Individual 0.17 [0.04, 0.31]
Year 0.03 [0, 0.12]
Residual 0.42 [0.29, 0.55]
Adjusted repeatability 𝑟
Individual 0.28 [0.08, 0.46]
Year 0.04 [0, 0.17]
Residual 0.68 [0.49, 0.89]
Sample sizes 𝑛
Years 3
Individuals 83
Observations 164
Code
plot(allEffects(mod_max_score))

age- and sex-specific variation in breeding plumage (of known-aged birds)

Code
# draw gt table
mod_max_score_age_sex_table <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = c(estimate),
             rows = 1:12,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = c(estimate),
             rows = 11:15,
             decimals = 0,
             use_seps = FALSE) %>% 
  sub_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = c(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

mod_max_score_age_sex_table
Banded Dotterel rufous band score Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Male (Sex) 1.09 [0.59, 1.57]
Age 0.04 [-0.21, 0.26]
Partitioned 𝑹²
Total Marginal 𝑹² 0.40 [0.19, 0.64]
Sex 0.35 [0.12, 0.6]
Age 0.00 [0, 0.3]
Total Conditional 𝑹² 0.66 [0.41, 0.86]
Random effects 𝜎²
Individual 0.17 [0, 0.41]
Year 0.04 [0, 0.2]
Residual 0.26 [0.1, 0.43]
Adjusted repeatability 𝑟
Individual 0.35 [0, 0.7]
Year 0 [0, 0.33]
Residual 1 [0.23, 1]
Sample sizes 𝑛
Years 3
Individuals 24
Observations 42
Code
plot(allEffects(mod_max_score_age_sex))

migratory status- and sex-specific variation in breeding plumage (of known-status birds)

Code
# draw gt table
mod_max_score_table_status_sex <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = c(estimate),
             rows = 1:10,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = c(estimate),
             rows = 11:13,
             decimals = 0,
             use_seps = FALSE) %>% 
  sub_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = c(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

mod_max_score_table_status_sex
Banded Dotterel rufous band score Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Male (Sex) 1.23 [0.82, 1.63]
Resident (Migratory Status) 0.48 [0, 0.92]
Partitioned 𝑹²
Total Marginal 𝑹² 0.51 [0.31, 0.69]
Sex 0.50 [0.31, 0.69]
Migratory Status 0.06 [0, 0.34]
Total Conditional 𝑹² 0.59 [0.38, 0.79]
Random effects 𝜎²
Individual 0.06 [0, 0.25]
Residual 0.30 [0.15, 0.45]
Adjusted repeatability 𝑟
Individual 0.16 [0, 0.54]
Residual 0.84 [0.46, 1]
Sample sizes 𝑛
Years 3
Individuals 21
Observations 42
Code
plot(allEffects(mod_max_score_status_sex))

within-individual moult dynamics (sex)

Code
# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in sex while accounting for sex differences in maximum rufous band score.
ind_prop_molt_scores <- 
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score))

# Assess sample sizes of each sex
ind_prop_molt_scores %>% 
  group_by(sex) %>% 
  summarise(n_distinct(rings_comb))
# A tibble: 2 × 2
  sex   `n_distinct(rings_comb)`
  <chr>                    <int>
1 F                           49
2 M                           35
Code
# Step 1 and 2: Standardize max_breeding_score within each sex group
ind_prop_molt_scores <- 
  ind_prop_molt_scores %>%
  group_by(sex) %>%
  mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>%
  ungroup()

# mixed effects binomial model comparing sex and date effect on the changes in moult scores
mod1 <- 
  lme4::glmer(prop_molt_score ~ 
                date_J * sex + std_max_breeding_score +
         # date_J * sex +
         (1 | rings_comb) + (1|season),
       data = ind_prop_molt_scores, 
       family = binomial,
       control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

plot(allEffects(mod1))

Code
# strong date effect, but no sex effect
tbl_regression(mod1, intercept = TRUE, 
               label = list(date_J ~ "Date", sex ~ "Sex"))
Characteristic log(OR)1 95% CI1 p-value
(Intercept) 20 15, 24 <0.001
Date -0.09 -0.11, -0.07 <0.001
Sex
    F
    M 1.1 -5.7, 8.0 0.8
std_max_breeding_score -0.29 -0.65, 0.07 0.11
Date * Sex
    Date * M -0.01 -0.04, 0.03 0.7
1 OR = Odds Ratio, CI = Confidence Interval
Code
# extract predicted trends
pr <- ggeffects::predict_response(mod1, c("date_J [30:293]", "sex"))

# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(date_J = x,
         sex = group) %>% 
  left_join(., dates_for_plot, by = "date_J")

# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = date, y = predicted, color = sex)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex),
              lwd = 1, alpha = 0.25) +
  luke_theme +
  theme(legend.position = c(0.3, 0.2),
        legend.justification = c(1, 0),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 10, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1)) +
  ylab("proportion intact of individual-based maximum rufous band") +
  scale_x_date(date_labels = "%B", 
               expand = c(0.01, 0.01), 
               date_breaks = "1 month") +
  scale_colour_brewer(palette = "Dark2", direction = -1,
                      name = "Sex",
                      labels = c("Female (N = 49)", "Male (N = 35)")) +
  scale_fill_brewer(palette = "Dark2", direction = -1,
                    name = "Sex",
                    labels = c("Female (N = 49)", "Male (N = 35)"))

within-individual moult dynamics (migratory status)

Code
# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.
ind_prop_molt_scores <- 
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, migratory_status, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score)) %>% 
  filter(migratory_status %in% c("R", "M"))

# Assess sample sizes of each sex
ind_prop_molt_scores %>% 
  group_by(sex, migratory_status) %>% 
  summarise(n_distinct(rings_comb))
# A tibble: 4 × 3
# Groups:   sex [2]
  sex   migratory_status `n_distinct(rings_comb)`
  <chr> <chr>                               <int>
1 F     M                                       2
2 F     R                                      11
3 M     M                                       3
4 M     R                                       5
Code
# mixed effects binomial model comparing sex and date effect on the changes in moult scores
mod1_status <- 
  lme4::glmer(prop_molt_score ~ 
                date_J * migratory_status +
         # date_J * sex +
         (1 | rings_comb) + (1 | season),
       data = ind_prop_molt_scores, 
       family = binomial,
       control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

# strong date effect, but no status effect
tbl_regression(mod1_status, intercept = TRUE, 
               label = list(date_J ~ "Date", migratory_status ~ "Migratory status"))
Characteristic log(OR)1 95% CI1 p-value
(Intercept) 25 10, 40 0.001
Date -0.11 -0.19, -0.04 0.002
Migratory status
    M
    R 12 -17, 41 0.4
Date * Migratory status
    Date * R -0.06 -0.20, 0.08 0.4
1 OR = Odds Ratio, CI = Confidence Interval
Code
# extract predicted trends
pr <- ggeffects::predict_response(mod1_status, c("date_J [30:293]", "migratory_status"))

# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(date_J = x,
         migratory_status = group) %>% 
  left_join(., dates_for_plot, by = "date_J")

# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = date, y = predicted, color = migratory_status)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status),
              lwd = 1, alpha = 0.25) +
  luke_theme +
  theme(legend.position = c(0.3, 0.2),
        legend.justification = c(1, 0),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 10, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1)) +
  ylab("proportion intact of individual-based maximum rufous band") +
  scale_x_date(date_labels = "%B", 
               expand = c(0.01, 0.01), 
               date_breaks = "1 month") +
  scale_colour_brewer(palette = "Set1", direction = -1,
                      name = "Migratory status",
                      labels = c("Resident (N = 16)", "Migratory (N = 5)")) +
  scale_fill_brewer(palette = "Set1", direction = -1,
                    name = "Migratory status",
                      labels = c("Resident (N = 16)", "Migratory (N = 5)"))

within-individual moult dynamics (age)

Code
# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.
ind_prop_molt_scores <- 
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, age, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score)) %>% 
  filter(!is.na(age))

# Assess sample sizes of each sex
ind_prop_molt_scores %>% 
  group_by(age) %>% 
  summarise(n_distinct(rings_comb))
# A tibble: 6 × 2
    age `n_distinct(rings_comb)`
  <dbl>                    <int>
1     1                        6
2     2                        5
3     3                       12
4     4                       13
5     5                        4
6     6                        2
Code
# mixed effects binomial model comparing sex and date effect on the changes in moult scores
mod1_age <- 
  lme4::glmer(prop_molt_score ~ 
                date_J + age +
         # date_J * sex +
         (1 | rings_comb),# + (1 | season),
       data = ind_prop_molt_scores, 
       family = binomial)
       # control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

# strong date effect, but no status effect
tbl_regression(mod1_age, intercept = TRUE, 
               label = list(date_J ~ "Date", age ~ "Age"))
Characteristic log(OR)1 95% CI1 p-value
(Intercept) 53 18, 88 0.003
Date -0.28 -0.47, -0.10 0.003
Age 2.2 0.61, 3.7 0.006
1 OR = Odds Ratio, CI = Confidence Interval
Code
# extract predicted trends
pr <- ggeffects::predict_response(mod1_age, c("date_J [30:293]", "age"))

# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(date_J = x,
         age = group) %>% 
  left_join(., dates_for_plot, by = "date_J")

# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = date, y = predicted, color = age)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = date, ymax = conf.high, ymin = conf.low, fill = age),
              lwd = 1, alpha = 0.25) +
  luke_theme +
  theme(legend.position = c(0.3, 0.2),
        legend.justification = c(1, 0),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 10, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1)) +
  ylab("proportion intact of individual-based maximum rufous band") +
  scale_x_date(date_labels = "%B", 
               expand = c(0.01, 0.01), 
               date_breaks = "1 month", 
               limits = c(as.Date("2021-10-01"), as.Date("2022-05-01"))) +
  scale_colour_brewer(palette = "Spectral", direction = -1,
                      name = "Age",
                      labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", 
                                 "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) +
  scale_fill_brewer(palette = "Spectral", direction = -1,
                    name = "Age",
                      labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", 
                                 "4 (N = 13", "5 (N = 4)", "6 (N = 2)"))

13-Nov-24: meeting notes (Bart, Luke, and Bashar via zoom)

  • add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.

  • drop pre-October data from the molt timing analysis.

  • next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.